Manifestation of Resonance-Related Chaos in Coupled Josephson Junctions 
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Chaotic features of systems of coupled Josephson junctions are studied. Manifestation of chaos in 
the temporal dependence of the electric charge, related to a parametric resonance, is demonstrated 
through the calculation of the maximal Lyapunov exponent, phase-charge and charge-charge Lis- 
sajous diagrams and correlation functions. The number of junctions in the stack strongly influences 
the fine structure in the current voltage characteristics and a strong proximity effect results from 
the nonperiodic boundary conditions. The observed resonance-related chaos exhibits intermittency 
over a range of conditions and parameters. General features of the system are analyzed by means of 
a linearized equation and the criteria for a breakpoint region with no chaos are obtained. Such cri- 
teria could clarify recent experimental observations of variations in the power output from intrinsic 
Josephson junctions in high temperature superconductors. 

PACS numbers: 74.81. Fa, 74.50.+r, 74.40.De, 74.78.Fk 



I. INTRODUCTION 

Systems of coupled Josephson junctions (JJs) are 
prospective objects for superconducting electronics 
and therefore continue to be the subject of intense 
investigations* 1 - - — Intrinsic Josephson effects - tunneling 
of Cooper pairs between superconducting layers inside 
of strongly anisotropic layered high temperature super- 
conductors (HTSCs) - provide an opportunity to model 
HTSCs as systems of coupled intrinsic Josephson junc- 
tions (IJJs) and an effective way to investigate nonlinear 
effects and non-equilibrium phenomena in HTSCs. In- 
trinsic tunneling is a powerful method for the study of 
the nature of high temperature superconductivity, trans- 
port along a stack of superconducting layers, and the 
physics of vortices. It also plays an important role in de- 
termining the current-voltage characteristics (CVCs) of 
tunneling structures based on HTSCs and the properties 
of vortex structures in these materials. 

Parametric resonance in coupled JJs demonstrate a se- 
ries of novel effects predicted by numerical simulations. 10 
In comparison with the single junction, the system of cou- 
pled JJs has a multiple branch structure in its CVCs. The 
outermost branch of the CVCs has a breakpoint (BP) 
and a breakpoint region (BPR) before a transition to an- 
other branch, caused by parametric resonance.— The BP 
current characterizes the resonance point at which a lon- 
gitudinal plasma wave (LPW) is created. It was demon- 
strated that the CVC of the stack exhibits a fine structure 
in the BPRi 1 ^ The breakpoint features were predicted 
theoretically^ 3 - and recently observed experimentally— A 
breakpoint region also appears in numerical simulations 
performed by other authors. 15 

The physics of IJJs, which are nonlinear systems, can- 
not be completely understood without the investigation 
of their chaotic features. Parametric instabilities are 



common phenomena in one-dimensional arrays of Joseph- 
son junctions. For example, a one-dimensional parallel 
array of N identical Josephson junctions was studied by 
discrete sine-Gordon equation (also known as Frenkel- 
Kontorova model) in Reffl6l. The parametric instabili- 
ties were predicted by theoretical analysis and observed 
experimentally. In particular, the novel resonant steps 
related to the parametric instability were found in the 
current voltage characteristics of discrete Josephson ring. 
It was verified experimentally that such steps occur even 
if there were no vortices in the ring. A long Josephson 
junction is also often considered as a one-dimensional 
parallel array of N identical Josephson junctions. An 
experimental testing of the parametric instability in IJJ 
would clarify the generic feature of one-dimensional par- 
allel and series array models and the role played by res- 
onance related chaos in these systems. 

The chaotic features of IJJs in HTSCs are also of great 
interest due to the observed powerful coherent THz radi- 
ation from such systems^ Broadly tunable sub-terahertz 
emission was found in the CVC near the switching point 
from the outermost branch to inner branches, and also 
from inner branches £ We stress that this is the same 
region in the CVCs where the BP and BPR were ob- 
served. Recently, for example, coherent THz radiation 
has been measured experimentally and was associated 
with a "bump" structure corresponding to the same part 
of the CVCsj 17 i 18 Obviously, the chaos in IJJs could 
strongly effect their radiation properties, because the 
transition to the chaotic state may inhibit the coherent 
radiation from IJJs. Therefore, detailed investigations of 
such transitions in real systems are very important. This 
situation makes the phase dynamics and the investiga- 
tion of chaos in intrinsic JJs, near the BPR of the CVCs, 
an urgent problem. 

In an early study of the rf current driven Josephson 
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junction it was found that, depending on the amplitude 
of the rf current, the system would develop a chaotic 
character.— ~— It was also shown that the onset of chaos 
as a function of frequency correlated with the function 
that indicated infinite gain, also as a function of fre- 
quency, for the unbiased parametric amplifier.— More 
recently the importance of chaos in IJJs, and its effects 
on the CVCs and the Shapiro steps in these systems, was 
stressed in Refs. HUH 

In a further comprehensive study, the critical behavior 
of the dynamical equation, in the sense of how it goes 
from regular to chaotic, was investigated^ It was dis- 
covered that the critical behavior is that of the circle 
map. Here, in the subcritical state, the resonances are 
separated by quasi-periodic orbits (i.e. having irrational 
winding number), whereas in the supercritical state the 
dynamics is composed of chaotic jumps between reso- 
nances. The chaos appears to be the result of the overlap 
of resonances.— 

In a similar study, the interaction of the Josephson 
junction with its surroundings was considered by cou- 
pling it in parallel with an RLC circuit, modeling a res- 
onant cavity^ The chaotic nature develops through the 
familiar period doubling sequence of bifurcations. Other 
studies exist that try to understand the irregular response 
of dissipative systems driven by external sources when 
such behavior is interleaved by the synchronized motion 
(Shapiro steps). This so called intermittent chaos is mod- 
eled by a random walk between the two neighboring steps 
that have become unstable. It is shown that the power 
spectrum of the voltage correlations has a broadband 
background characterizing the chaotic solution^ It may 
be noted that a similar random walk model has been de- 
veloped to explain the phase synchronization of chaotic 
rotators^ 

Interesting features of chaos, chaos control and chaotic 
synchronization in Josephson junction arrays ( J JAs) and 
shunted JJ models have been described in Refs. l3ll - 
JJAs attract much attention due to their perspec- 
tive as high precision voltage standards and, as men- 
tioned earlier, high power coherent THz sources. Out- 
put power of a single junction is extremely low (typically 
10 nW), but much higher output power can be obtained 
by using JJAs. Chaos and nonlinearity in JJAs have 
been described in Refs. l34l|35l In particular, Basler et 
al£5. reported on the theory of phase locking and self- 
synchronization in a small JJ cell by using the resistively 
shunted junction (RS J) model and Hassel et al£^ studied 
self-synchronization in distributed Josephson junction ar- 
rays. Features of spatiotemporal chaos in JJAs composed 
of RS Js were numerically investigated by Bhagavatula et 
aL— Chernikov and Schmidt 38 demonstrated that when 
four or more JJs are globally coupled by a common re- 
sistance, the system exhibits adiabatic chaos. Direct cal- 
culation results 34 showed that phase locking and chaos 
coexist when there are three junctions in the JJA. In 
Ref. [33| chaos, hyperchaos and controlling hyperchaos in 
a JJA were investigated. The numerical results show that 



chaos and hyperchaos states can coexist in this array of 
two shunted JJs, independently of whether or not the 
original states of the junctions are chaotic. 

Although there have been recent reports on the 
chaotic behavior in JJs, such as Ref. |3J, there are 
relatively few detailed investigations of chaotic behav- 
ior in closely-related phenomenological models such as, 
the capacitively-coupled model (CCJJ),— the resistive- 
capacitive shunted model (RCSJJ)) 31 ' 40 or the capaci- 
tively coupled Josephson junction model with diffusion 
current (CCJJ+DC) 41 ' 42 of the present work. Specifi- 
cally, the chaotic features at the parametric resonance 
have not been investigated before. 

In this paper, we study the phase dynamics in coupled 
Josephson junctions. A resonance-type chaos related to 
the parametric resonance in this system is demonstrated. 
The origin of chaos is related to the coupling between 
junctions based on the fact that the superconducting lay- 
ers (S-layers) are in the nonstationary nonequilibrium 
state due to the injection of quasiparticle and Cooper 
pairs i 39 i 43 We present results on a chaotic mode which 
is due to the coupling between junctions and a paramet- 
ric resonance in coupled JJs. This mode cannot be ob- 
served in case of a single J J. Creation of a LPW and 
its interplay with a discrete lattice of superconducting 
layers produces a rich dynamic behavior, including, for 
example, intermittency. We use different methods to in- 
vestigate this dynamics, including correlation functions 
and Lyapunov exponents. Results for stacks with differ- 
ent numbers of JJs, in which LPWs with different wave 
numbers are exited, give interesting information on cer- 
tain features of the fine structure within the BPR of the 
CVCs. The influence of the boundary conditions on the 
chaotic part of the fine structure leads to the manifes- 
tation of a strong proximity effect in IJJs. Transitions 
between chaotic and regular behavior in BPR are pre- 
dicted. Transitions from the chaotic state to both a reg- 
ular LPW state and a traveling wave mode are demon- 
strated. Finally some general features of the system are 
analyzed and the conditions and parameters leading to a 
BPR with no chaos are found. 

The paper is organized in the following way. In Sec- 
tion II we introduce the coupled sine-Gordon equation, 
used for numerical simulations, and discuss the bound- 
ary conditions and numerical procedures. We briefly de- 
scribe the method for calculating the Lyapunov expo- 
nent. Section III is devoted to the study of the resonance- 
related chaos produced in the temporal dependence of 
the electric charge, and this is demonstrated through the 
calculation of the maximal Lyapunov exponent, phase- 
charge and charge-charge Lissajous diagrams and corre- 
lation functions. We also provide examples of the man- 
ifestation of related chaotic features in the correlation 
functions and polar diagrams. The effect of the num- 
ber of junctions in the stack is considered in Sec. IV. 
Nonperiodic boundary conditions and the proximity ef- 
fect in coupled JJs are discussed in Sec. V. We show the 
occurrence of transitions from chaotic to regular behav- 
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ior in Sec. VI. To clarify the chaotic features in coupled 
JJs, we analyze the linearized equation for phase differ- 
ences in Sec. VII, and discuss the question concerning the 
region of instability at chaotic boundary. In Sec. VIII 
we demonstrate the parametric resonance region without 
chaotic behavior. Conclusions follow in Sec. IX. 



II. MODEL AND METHOD 

High-Tc superconductors have a layered crystal struc- 
ture and strong anisotropy. A model of multilayered 
Josephson junctions, which has an alternating stack of 
superconducting and insulating layers, is sufficient to de- 
scribe the electronic properties of these systems^. In the 
present model we assume that the physical quantities are 
spatially homogeneous on each layer. This assumption is 
applicable in cases of no applied magnetic field and for a 
small sample size in the 06-direction, i.e. for the inten- 
sively developing research field known as the physics of 
Josephson nano junctions, where the length of the Joseph- 
son junctions is smaller than the Josephson penetration 
length. One-dimensional models with coupling between 
junctions do capture the main features of real I JJs, like 
hysteresis and branching of the current voltage charac- 
teristics, and thus help to improve the understanding of 
their physics. An interesting and very important fact is 
that the ID models can also be used to describe the prop- 
erties of a parallel array of Josephson junctions, which is 
often considered as a model for long Josephson junctions. 
The experiments and theoretical studies on the propaga- 
tion of Josephson fluxons and electromagnetic waves in 
parallel arrays of Josephson junctions were presented in 
Refs.? ? . Their experimental data demonstrate a series 
of resonances in the current voltage characteristics of the 
array and were analyzed using the discrete sine- Gordon 
model and an extension of this model which includes a ca- 
pacitive interaction between neighboring Josephson junc- 
tions. 

To simulate the CVCs of intrinsic Josephson junctions 
in high temperature superconductors we investigate the 
phase dynamics within the framework of the CCJJ+DC 
model.— ~— In this model the dynamical equations are 
given by, 



dV, 
dt 



=V t -a(V l+1 



I + I n - sin ipi - 
Vi-i 



dt ' 



(la) 
(lb) 



Here Vi = Vi+u is the voltage difference between 
the (I + l)th and Zth S-layers, and the gauge-invariant 
phase difference <pi between adjacent S-layers is given by 

= 0,+i(t) - 0,(t) - 1 SHtl )d dzA z {z,t), where 9 l 
is the phase of the order parameter (i.e. phase of the 
macroscopic wave function of the superconductor) in the 
Ztb. S-layer of thickness s, A z is the vector potential in 
the barrier of thickness D, and d = s + D is the period of 



the lattice. In the above equations, a and (3 are the cou- 
pling constant and dissipation parameter, respectively. 
Time t is normalized to the inverse plasma frequency 
uip 1 (with ujp = 2eI c /HC, and C denoting the junction 
capacitance) , the voltage differences Vi are normalized to 
Vo = ftio p /(2e), and the bias current I and added noise 
current /„ are normalized to the critical current I c . The 
role of the added noise current (/„ ~ 10 -8 ) has been dis- 
cussed previously^ We solve this system of dynamical 
equations for stacks with various numbers (N) of JJs. 
The system of Eqs. ([T]) can also be written in the form 



dt 2 



= £ ^ 



I + I n — sin ip v - (3 



dipt 



dt 



(2) 



where, for nonperiodic boundary conditions (BCs), the 
matrix A has the form 



A = 



1 + aG —a 
~a l + 2a 
-a 1 





-a ... 

I- 2a —a 



V 



-a 1 



\ 



xGJ 



(3) 



Here I and V run over the N junctions, (3 C — lo 2 R 2 C 2 is 
the McCumber parameter, and G = 1 + 7, 7 = s/sq — 
s/sn, where s, sq and sn are the thicknesses of the mid- 
dle, first, and last S layers, respectively. According to 
the proximity effect, the thicknesses of the first and last 
layers are assumed to be thicker than the middle layers 
inside the stack. So, for the nonperiodic boundary BCs, 
the equations for the first and last layers in Eq. ^ are 
different from those of the middle layers i 39 i 46 For periodic 
BCs the matrix A has the form 



(\ + 2a 



A = 



—a 




-a 
- 2a —a 
-a l + 2a 



-a 



V 



(4) 



-a l + 2a/ 



Using the Maxwell equation eeodivE = Q, where e 
is the dielectric constant of the insulating layers and Eq 
is the permittivity of free space, we express the charge 
density in the Zth S-layer as 



Qi = Q a(Vi +1 - Vi), 



(5) 



where Qq = eeoVo/rfj and rjj is the Debye screening 
length. In all our calculations the CVCs and time de- 
pendence of the charge oscillations in the S-layers are 
simulated at a = 1. The system not very sensitive to the 
value of the coupling parameter. Further details concern- 
ing numerical simulations of the CVCs and the electric 
charge can be found in Refs. I46l - l48l 

The usual test for chaos in a system is through the cal- 
culation of the largest Lyapunov exponent (LE)4^ The 
general idea is to follow two nearby trajectories and to 
calculate their average logarithmic rate of separation. 
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Whenever the two trajectories get too far apart, one of 
them has to be moved back to the vicinity of the other, 
along the line of separation. A conservative procedure is 
to do this at each iteration. 

In this paper we solve the system of Eqs. (UJ, using 
a fourth-order Runge-Kutta method, to find the V\ and 
(pi. To this end we choose the time domain Tf and time 
step T p — Tf/m, such that tj = jT pi for j = 0, . . . ,m. 
In order to calculate the largest LE, at fixed current I, 
we first integrate the system for a time Ti, until all tran- 
sients have decayed. Then, starting from an arbitrary 
small initial separation vector dj = (AV(tj), A(p(tj)), we 
advance both systems by one time step and calculate 
the new separation dj+i = (AV(tj + i), Aip(tj + i)), where 
A<pi(U) = (ff^U) ~ <Pi(U) and AVi{U) = V{{U) - Vi(U). 
The Lyapunov exponent for this step is defined by 



LEi = — In 



(6) 



where \\d,\\ = E^A^fo)) 2 + (A^fe))T /2 - In or- 
der to avoid cumulative numerical round-off errors, due 
to the expected exponential changes in separation, the 
separation vector is rescaled after each step in time, i.e. 
we set dj + i = (||dj||/||d.,- + i||)(ij + i, before advancing to 
the next step. This procedure is iterated over the rest 
of the time domain at fixed current. Finally we take the 
average of all the LEj to find the Lyapunov exponent 



LE(I) 



1 



m— 1 

■E 

j=k 



LE, 



(7) 



where tk = kT p = Ti. By changing the current in small 
steps (d/), and repeating the above procedure for each 
value of current, the LE is obtained as a function of 
current. Typically, except where noted, we have used 
Ti = 50, T f = 25000, m = 5 x 10 5 , and dl<5x 10~ 4 . 
Additional details concerning the LE calculation can be 
found in Refs. SiHH 



First, we consider a stack of 9 JJs using periodic BCs 
with /3 = 0.2. In Fig. [TJ the Lyapunov exponent and 
outermost branch of CVCs, as functions of current, are 
compared with the charge oscillation in the 8th S-layer, 
as a function of time and current. 

To calculate the voltages Vi(I) at each point of CVC 
(i.e. for each value of /), we simulate the dynamics of the 
phases <pi (t) by solving the system of equations (TJJ using 
the fourth-order Runge-Kutta method with a time step 
T p . After simulation of the phase dynamics we calculate 
the dc voltages on each junction by averaging over a cho- 
sen time domain. After completing the voltage averaging 
for the fixed bias current /, the bias current is increased 
or decreased by a small amount dl (bias current step) 
in order to calculate the voltages in all junctions at the 
next point of the CVC. We use the distribution of phases 
and their derivatives at the end of the time domain of 
previous point as the initial distribution for calculating 
the next point. In this manner the dependence of charge 
on current 'includes' the time dependence of each current 
step. 

In the inset we illustrate the position of the BPR on 
the outermost branch of the CVCs. As one can observe, 
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III. INDICATIONS OF CHAOS IN THE BPR 
FOR A SYSTEM OF COUPLED JJs 

The CVCs of coupled JJs has a multiple branch struc- 
ture. Here we concentrate on the outermost branch only. 
The system of coupled JJs exhibits a parametric reso- 
nance which is reflected in the CVC as a breakpoint. In 
stacks with odd numbers of junctions a breakpoint re- 
gion, with a fine structure corresponding to the width 
of the parametric resonance, is observed. A study of the 
fine structure in the BPR gives one enough reason to look 
for chaos in the dynamics of a stack. Its manifestation 
depends strongly on the parameters of the system. Here 
we present some examples and methods, which might be 
used for demonstration and investigation of the chaotic 
features of coupled JJs. 



FIG. 1: (Color online) Lyapunov exponent and CVC, as func- 
tions of current, with the charge oscillation in Q% (the charge 
in the 8th S-layer) as a function of time. Inset demonstrates 
the position of the BP on the outermost branch of the CVCs. 

there are two distinct regions: LE — and LE > 0. The 
case LE > strictly indicates that the system is chaotic 
over the approximate current interval of (0.5620, 0.5536), 
while LE = indicates marginal stability for the system 
elsewhere. Both the CVCs and charge oscillations show 
chaotic features which coincide with LE > 0. 

Second, we study the correlations in superconduct- 
ing currents in the neighboring Josephson junctions, and 
charge correlations in neighboring superconducting lay- 
ers. This allows us to further demonstrate the observed 
chaotic features in CVC. The charge-charge correlation 
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functions for neighboring layers are given by 

Ct l+1 {I) = {Qi(t)Qi +1 (t)) (8) 

= lim — l — [ f Qi(t)Qi+i{t)dt, 

7/-Kx> If — ±i J T . 

where I is the layer number and Ti is the initial time for 
the averaging procedure. Figure [5] presents the results of 
calculations for the cases I = 2 and I — 7, for the same 
stack of nine JJs. The negative sign of Cf l+1 is due to 
the fact that a LPW with wave number k = 87r/(9d) is 
created. This wave number is close to the 7r-mode, for 
which the charge on neighboring layers differs in sign, 
and hence the product of positive and negative charges 
produces a negative sign (cf. Sec. IV). Here the expected 
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FIG. 2: (Color online) The charge correlation functions 
Cm+i = (Ql{t)Ql+i(t)) for I = 2 (labeled by 2,3) and I = 7 
(labeled 7,8) for the stack with nine IJJ as a function of the 
bias current. CVC is also shown (and labeled by an arrow). 
A chaotic part of the BPR (encircled) is enlarged in the inset. 

feature of the correlation functions in the chaotic region 
is confirmed: at the transition to chaotic behavior (point 
C\), the values of all correlation functions approach each 
other, i.e. (Qi(t)Qi +1 {t)) = (Qi(t)) (Qi +1 (t)) . Thus, 
within the chaotic region, all correlations are lost. To 
demonstrate the character of the correlation function, we 
enlarged a region in the chaotic part of the BPR in the 
inset to Fig. [2] We see that both correlation functions 
6*2,3 and C73 show a similar variation with current, in 
agreement with the Lyapunov exponent. 

Third, manifestations of the chaotic behavior can also 
be demonstrated by phase-charge diagrams. To do so, we 
record the time dependence of the phase difference across 
the chosen junction and electric charge in the chosen S- 
layer at a given bias current. Figure [3] shows such a 
phase-charge diagram: the variation of the absolute value 
of charge in the first superconducting layer Qi with the 
phase difference <p\ across the first JJ of the stack. This 




(a) : 1=0.57126; 

(b) : 1=0.560; 



FIG. 3: (Color online) The phase-charge diagram. It shows 
the absolute value of the charge on the first superconducting 
layer (radial coordinate), versus the phase difference across 
the first junction at two values of bias current: (a) I — 0.57126 
and (b) I = 0.56000. 

is a Lissajous figure, showing that the motion at I = 
0.57126 is qualitatively different from that at I = 0.560, 
thus providing additional evidence of chaos in this part 
of the CVCs. 



IV. EFFECT OF THE NUMBER OF 
JUNCTIONS IN THE STACK 

We next study the behavior of the stacks as a function 
of the number of junctions. To clarify this effect, we first 
discuss the parametric resonance features in the stacks 
with even and odd numbers of JJs. 

In the stacks with an even number of junctions and 
periodic BCs, at j3 = 0.2, a LPW with wave number 
k = 7r (7r-mode) is created, so that its wavelength A = 2d 
satisfies the commensurability condition L = nXM- Here 
L = Nd is the length of the stack, and d is the period of 
the lattice. In this case a fine structure in the BPR is ab- 
sent and the transition from the outermost branch to an- 
other one is observed after a sharp (exponential) increase 
of charge in S-layers. Also, the value of the breakpoint 
current is the same for all stacks with an even number of 
junctions. The sharp increase of charge in the S-layers 
happens in a very narrow interval of the bias current. 
For example, in a simulation of the CVC in a stack with 
N = 10 and a current step of 10 -4 , the interval is (0.5771, 
0.5767) in which the increase of charge is of about 7 or- 
ders of magnitude. 

In the stacks with an odd number of junctions the wave 
number of the LPW is k = (N — l)n/(Nd) which ap- 
proaches k = ir/d as N increases. The value of the BP 
current is limited by its value at k — 7T, and the width 
of the BPR decreases with increasing NM- With an odd 
number of junctions in the stack, the sharp increase of 
the charge accompanies the BP but then, before transi- 
tion to another branch, a wide BPR exists^ The origin 
of it is explained by the fact that the wavelength of LPW 
does not satisfy the commensurability condition. The 
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fine structure is related to the character of the charge 
oscillations in the superconducting layers^. A part of 
BPR reflects the chaotic behavior of the system, which 
is the main subject of our paper. Next we look at how 
N affects the structure of BPR. As mentioned above, 
the study of the correlation of charge on neighboring S- 
layers is a powerful tool to investigate the dynamics of 
charge in coupled JJs, within the BPR^ Here we use it 
to demonstrate the effect of the number of junctions on 
the structure of the BPR. 

For the periodic BCs the number assigned to a junc- 
tion is merely a label. However, in order to compare the 
correlation functions between stacks with different num- 
bers of junctions, in different simulations, it is useful to 
use a specific layer as reference. In the parametric reso- 
nance region the phase dynamics in the coupled system 
of JJs singles out such specific layers, namely, if a node of 
the charge standing wave coincides with the layer, then 
the charge on that layer is close to zero. Therefore we 
number the correlation functions according to their dis- 
tance from such specific layers, as shown in Fig. 2) To 



9 

8 ~ 

7 ~ 

6 ~ 

5 ~ 

4 Specific layer 

3 ~ 

2 ~ 

1 ~ 
9 



C 

C 4+ 



5+ 
c 



FIG. 4: (Color online) Formation of the charge correlation 
functions in a periodic stack of 9 IJJs. The numbers count 
the layers in the stack. 

reflect the symmetry of the correlation function relative 
to the chosen specific layer, we introduce the subscripts 
+ and — on the correlation functions, as shown in Fig. [4] 
In this case the specific layer is number 4, and we have 



C3 4 and so on. 



c i+ = c h and c t 

Figure [5] shows the part of the outermost branch of 
CVC in the BPR together with two correlation functions 
C% + and Cf_ for the stacks with 7, 9, 11 and 13 JJs with 
periodic BCs, at /3 = 0.2. In each case the part around 
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FIG. 5: (Color online) The CVC of the outermost branch 
in BPR together with correlation functions C%+ and C%- for 
periodic stacks with different numbers of IJJs: (a) N — 7; (b) 
N = 9; (c) N = 11; (d) TV = 13. 



B c corresponds to the maximal increase of charge in the 
S-layers (c.f. the beginning of charge dependence in Fig. 
2) and is characterized by a LPW frequency ujlpw- To 
distinguish different parts of BPR, we introduce the no- 
tation: B for the beginning of the parametric resonance, 
B c for the BP of the CVC, T and W to indicate the ap- 
pearance of low frequency charge modulations, and C\ to 
separate the chaotic part of the BPR. The loops that are 
seen in the correlation functions to the left of B c (very 
clearly manifested for N = 13), correspond to the charge 
dynamics on the specific layer in the regular part of the 
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BPR and were discussed in Ref. 

In the part T — Y a modulation of the charge oscil- 
lations with low frequency appears^ The part Y — C\ 
is the transition region between regular and chaotic be- 
havior where a second low frequency is observed. The 
features in the transition region Y — C\ relate to the spe- 
cific behavior of the LPW with k = (N—1)tt/(N<I) in the 
stacks with different N. We found that the region after 
C\ demonstrates the chaotic behavior of coupled system 
of Josephson junctions. 

From Fig. [5] one can draw the following conclusion: 
The increase in N increases the value of the breakpoint 
current, shifting the position of the BPR and decreasing 
the width of the BPR. As one can see, the increase in N 
retains the main features of the BPR structure for the 
stacks with 7, 9, 11 and 13 junctions. Simulation of the 
time dependence of the charge in S-layers in these regions 
confirms the same character of the charge oscillations for 
all these stacks. The detailed features seen in the fine 
structure of the correlation functions, for the stacks with 
different numbers of junctions, and hence also different 
lengths, are related to the specific behavior of LPWs with 
different wave numbers. As will be discussed in Sec. I VII 
transitions from chaotic to regular and back can be ob- 
served clearly in, for example, Fig. [5fa), at I ~ 0.556. 

Although the width of the BPR decreases with increas- 
ing N , the width corresponding to the chaotic behavior 
increases with increasing N . This fact is demonstrated 
in Fig. [6l where the LE as a function of / is presented for 
stacks consisting of different numbers of junctions. We 
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FIG. 6: (Color online) Maximal Lyapunov exponents for peri- 
odic stacks with odd numbers of junctions N. Double arrows 
emphasize the increase in width of the chaotic region, relative 
to the BPR, as N changes from 5 to 15. 

see that the chaotic part of the CVC with five junctions 
is significantly smaller that for 15 junctions. This fact 
has been stressed by the double arrows shown in Fig. [BJ 



V. PROXIMITY EFFECT 

Contact of the stack of superconducting layers with a 
normal metal (electrodes) leads to the proximity effect: 
the superconducting regions are expected to penetrate 
into the electrodes. Due to the proximity the effective 
thicknesses of the superconducting layers at the edges of 
the stack, sq and sn, are larger than the thicknesses, s, of 
the middle layers. As mentioned in connection with Eq. 
@, this effect can be taken into account by variation of 
the parameter 7 in the nonperiodic BCs4£ Therefore, in 
this section, we consider briefly the influence of 7 on the 
chaotic behavior of the CVC and correlation functions. 

The CVC and correlation function, Cf = Cf 6 (I) = 
(QsityQeft)) (I), at 7 = 0,0.5 and 1 are presented in 
Fig. [7] for a stack of 9 junctions. In this case, the LPW 
with k = (N — l)ir/N is created at the breakpoint for all 
stacks, with an even and odd number of junctions and 
with j3 = 0.2. 48 Figure [TJa) shows the case 7 = 0. As 
we see, the main features of the CVC and the correla- 
tion functions are coincident. Within the chaotic region 
there appears to be some regions with regular oscilla- 
tions, marked by arrows. Figure EJb) clearly shows two 
such regions in CVC and in the correlation function at 
7 = 0.5. With increasing 7 the width of the BPR is 
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FIG. 7: (Color online) The CVC and correlation function C c s 
in the BPR for the stack with 9 IJJ at (a) 7 = 0; (b) 7 = 0.5; 
(c) 7 = 1. Inset to (c) shows an enlarged view of the encircled 
area near the BP. 

decreased, especially the B c — T region, which becomes 
very small by the time 7 = 1. An enlarged view of this 



small region is shown in the inset to Fig. [Jjc). 

Taking into account the proximity effect further em- 
phasizes the influence of the number of junctions in the 
stack. In particular, for 7 = 1, the entire BPR becomes 
chaotic as N increases beyond about 15. This result is 
presented in Fig. where the CVC of the outermost 
branch together with LE are compared for N = 7 and 
N = 15. We see that in the case N = 15, the LE be- 



0.56 



0.565 



0.57 I 0.575 



18 
> 17 

16 
15 



N = 7, 
y=l 



CVC 



0.58 
- — I 0.06 



0.04 
0.02 




0.56 0.565 
0.57 

40 



0.57 1 0.575 
0.575 I 



0.58 
0.58 



39 



38 



37 ■ 





N= 15, 






y=l 

LE 


CVC 













0.06 
0.04 
0.02 




0.57 



0.575 



I 



0.58 



FIG. 8: (Color online) CVC of the outermost branch for stacks 
with (a) N = 7 and (b) TV = 15 JJs together with the LE at 
7 = 1. 

comes positive practically just after the breakpoint. In 
view of these results we predict that the proximity effect 
will have a strong influence on the fine structure of the 
BPR in the CVC of real junctions. 



VI. INTERMITTENCY IN CVC AND 
CORRELATION FUNCTIONS 

As we have seen in Fig. there are some regular re- 
gions within the chaotic parts of the CVC and correlation 
functions, and so transitions exist from chaotic to regu- 
lar behavior and back. This is known as intermittency 
in dynamical systems with chaotic dynamicsj 52 i 53 Now, 
for the first time, we demonstrate such intermittency for 
the resonance related chaos in systems of coupled JJs. 
Many such transitions can be seen in Fig. |H1 where we 
present results of a high precision calculation of the LE 
together with the CVC for a stack of 9 junctions, using 
nonperiodic BCs, at 7 = 0.5. In this calculation Ti = 50, 
T f = 250000, to = 5 x 10 6 and dl = 10~ 5 . The reg- 
ular regions seen in the CVC coinciding with LE = 0. 
In particular the two largest intervals with LE = 
correspond to the current intervals (0.5598, 0.5600) and 
(0.5573,0.5575), shown by double arrows. To better un- 
derstand the nature of these transitions, the charge time 
dependencies in these regions should be considered. In 




0.56 



0.57 



FIG. 9: (Color online) Intermittency is shown by windows 
of LE — inside the chaotic dynamics for a stack of 9 JJs 
with nonperiodic BCs. Arrows point to the respective scales 
for each curve and the dashed line shows the LE = axis. 
Double arrows show the correspondence between features of 
the CVC and the LE. 



future work it would also be of interest to calculate the 
full spectrum of Lyapunov exponents, which would allow 
one to distinguish between chaotic behavior (one pos- 
itive exponent) and hyperchaotic behavior (more than 
one positive exponent). 

To gain more information about the transitions we in- 
vestigate the dependence of all the charge-charge cor- 
relation functions and the LE, on the bias current in 
the BPR. As the LE shows in Fig. [TU1 the absence of 
the charge correlations in different S-layers is a signature 
of the chaotic behavior. Since we are interested in the 
chaotic features, we have not labeled each curve by its 
corresponding correlation function. In Fig. [TU] we see a 
restoration of correlations in the middle of the chaotic re- 
gion for a stack with 13 JJs. All presented characteristics 
(Cc, CVC and LE) reflect this transition from the chaotic 
behavior to regular and back. To stress the agreement 
between the correlation functions and LE, vertical dashed 
lines have been drawn in Fig. [TJJJ We see that changes in 
the LE (small peaks on the LE curve) correspond exactly 
to changes in the charge correlations. 

We have characterized one such transition from the 
chaotic to regular behavior as corresponding to a charge 
traveling wave stated As shown in Fig. [TTJ a transition 
from chaos to the traveling wave branch (TWB, dashed 
line) exists for the stack with 7 junctions, at /3 = 0.2. 
This transition is encircled in the figure. We can see in 
the upper inset that the time dependence of charge in 
this region demonstrates a changing of character of the 
charge oscillations in S-layer. In the other inset (bottom 
right of Fig. ITU we have shown the Lissajous charge- 
charge diagram in the chaotic state and in the TWB. 
The diagram demonstrate the variation in time of the 
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FIG. 10: (Color online) Demonstration of the intermittency 
within the interval of bias current 0.5659 < I < 0.5565 (shown 
by a double arrow) for the charge-charge correlation functions 
(Cc) and LE for a periodic stack of 13 junctions. 
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FIG. 11: (Color online) Transition in the outermost branch 
(OB) to the traveling wave branch (T WB) in the CVC of the 
stack with 7 JJ. Inset above the curve demonstrates charge vs. 
time for S-layer at the transition point. Inset below the curve 
shows the charge-charge diagram for the neighbor S-layers in 
the chaotic region of the OB and in TWB. 



charges in two neighboring layers. The charge in the first 
S-layer, Q\, is plotted along the x-axis, and the charge in 
the second S-layer, Q2, is plotted along the y-axis. As we 
can see, the Lissajous charge-charge diagram presents an 
open trajectory for the chaotic region and a closed one for 
the TWB. Chaos in systems of coupled JJs can also be 
transformed to regular behavior by external irradiation. 
This is effected by adding a time dependent current to I 
to the right hand side of Eq. (1), through the additional 
term, Asiri(ujRt), where A and lor are the amplitude and 
frequency of the microwave radiation. This result paves 



the way for experimental testing of the observed tran- 
sitions from chaotic to regular traveling wave behavior. 
Detailed investigations of the effects of radiation on IJ Js 
will be made in a future study. 

VII. LINEARIZED EQUATION: PARAMETRIC 
RESONANCE 

In this section we show that the linearized equation 
for the Fourier component Sk of the phase difference 81 = 
fi+i — (fi can determine the region of instability of the 
coupled JJ states, and that together with calculation of 
the Lyapunov exponent, in the framework of the same 
equation, we can arrive at some important conclusions 
concerning the dynamics of this system. In particular, 
there is a region of parameter values corresponding to the 
parametric resonance without transition to the chaotic 
state. This conclusion is supported by the results of the 
investigation of the temporal dependence of the electric 
charge in S-layers and the correlation functions. 

To investigate the instability region of coupled JJ 
states, corresponding to the outermost branch, we use the 
following linearized equation ! 11 ' 13 ? 48 It can be obtained 
from Eq. @ as, 

5\ + (l- aV (2) )(cos(^)<5/ + 05 t ) = 0, (9) 

by using the linear approximation sin^j+i) — sin(<p;) « 
Si cos(<fi), where <p ~ fit = Vt/N, il = V/N is the Joseph- 
son frequency, and V is the total voltage across the stack. 
Use of the Fourier expansion Si(t) — ^ fe <5/ c e 4 ' d , produces 

4 + P(k)5 k + cos(fL(k)t)5 k = 0, (10) 

where t is normalized to uj~ (k), with uj p (k) = 
w p C a , f3(k) = (3C a , Q{k) = fl/C a and C a = 
\Jl + 2a(l — cos(fc)). In Eq. (j9|), we have used the dis- 
crete Laplacian V (2) /i = fi+i + fi-i ~ 2/;. 

The linearized equation (|10[) demonstrates the insta- 
bility region on the /3(k) — Q(k) diagram, which char- 
acterizes the parametric resonance in the coupled JJj^ 
This linear equation is similar to a damped Mathieu 
type equation, with its periodic coefficient, demonstrat- 
ing parametric instability (resonance).— 

We calculate the LE for the dynamics dictated by 
the linearized equation (fTU|). by the method described 
above, for the same system considered in Fig. Q] Three- 
dimensional results are presented in Fig. Q21 In our cal- 
culations we have taken m — 300000, dfl(k) = 0.005 and 
df3(k) = 0.005. Below, we consider different cross sec- 
tions of this figure. The LE has positive values at small 
values of O(fc) and j3(k)) and is negative for large ft(k) 
and f3(k)). We also show the boundary between these 
regimes by the curve denoting LE = 0, in Fig [T3J 

To clarify the results presented in Fig. [121 we demon- 
strate two cross sections of it in the insets (a) and (b). 
Inset (a) shows the LE at constant value of dissipation 
parameter f3(k) = 0.7 as a function of f2(fe). We note that 
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(a) (b) 




FIG. 12: (Color online) Three dimensional picture of depen- 
dence of Lyapunov exponent versus fi(fc) and /3(k) given by 
equation (0. (a) The dependence of LE on Q(k) at fixed /3(k); 
(b) The dependence of LE on /3(k) at fixed fl(fc). 



the £l(k) is related to the voltage by Josephson relation 
= V I (NC(k)), i.e. f2(fc) is proportional to voltage. 
The dashed line shows LE = 0. For large enough il(fc), 
the LE is negative which points to regular dynamics in 
the system. For lower values of f2(fc) the behavior of the 
system is chaotic. Inset (b) shows the LE at constant 
value of voltage Sl(fc) = 0.7 as a function of /3(k). 

In Fig. [13] we show the cross section of Fig. [12] with 
LE = plane (solid line) together with the instability 
region (dashed line). Above the curve LE — 0, the be- 
havior of the system is regular, and below it, chaotic. The 
border of the instability region is related to the paramet- 
ric resonance in the system. The region between these 
two curves determines the regular behavior showing para- 
metric resonance. 

Comparison of the Fig. Q2] and Fig. Q] allows us to es- 
timate the width of the regular region at fixed f3(k), and 
check the accuracy of linearization and approximations 
we have used to derive Eq. (jlOp . From our definition of 
fl(k), we have AV = NCk£, where e = ASl(fc) denotes 
the width of regular region in direction of Q(k) axis (see 
Fig. HI?]) . In the case k = 8ir/9 (the corresponding BPR 
in CVC for the 9 coupled JJs with (3 = 0.2 is shown in 
Fig. □} we find C a = yjl + 2ev(l - cos(fc)) = 2.2. Then 
P(k) = (3C a = 0.44. From Fig. [13] for this (3(k) we have 
e = 0.071 and AV = NCkS ~ 1.4. As we can estimate 
from Fig. [T]the width of regular region is Ay « 1.6. So, 
the linearization and approximations used are accurate 




FIG. 13: (Color online) Lyapunov exponent (solid) and res- 
onance region (dashed) in Q(k)-/3(k) plane, e denotes the 
width of regular region in direction of fi(fc) axis. The hor- 
izontal double dashed dotted line corresponds to the value 
P{k) = 0.9772. A system with this value of j3{k) doesn't touch 
the chaotic region. Below we demonstrate that chaotic part 
of BPR does not appear for this value of /3(k) (see Fig. [KJ. 



enough and there is a good agreement between results 
presented in Figs. [L3l andfTI 



VIII. PARAMETRIC RESONANCE REGION 
WITHOUT CHAOTIC BEHAVIOR: CASE f3 = 0.6 

Let us now discuss in more detail the prediction which 
follows from Fig. [13] namely, that there is an interval of 
P{k) values for each fixed f2(fc) corresponding to the BPR 
without chaotic behavior. To find (3 which corresponds 
to a j3(k) in the discussed interval, we should know the 
wave number of the created LPW. We take an arbitrary 
value of /3(k) = 0.9772 in this interval. To determine the 
corresponding value of /3 we first find the wave number 
of LPW which is created in this case at the resonance 
point. For this purpose we investigate, in Fig. Q31 the /3- 
dependence of the breakpoint current I{, p (see details in 
Ref. fill ) in Fig. [H] Inset shows the enlarged part of the 
figure where changes of the wave number k of the LPW 
happen. We see that the value of (3 = 0.6 corresponds to 
k = 4tt/9. 

We examine our idea at f3 = 0.6 to check if for the 
stack with 9 coupled JJ a chaotic region is absent, and 
to see the character of charge oscillations in the S-layers. 
According to the relation (3(k) — f3C a , the value of j3 — 
0.6 follows from that chosen in Fig. [T3J namely /3(fc) = 
0.9772 and the LPW wave number k — Att/9 because 
C a = y/l + 2a(l - cos(fc)) = 1.6287. 

In Fig. [T5] we show the outermost branch in CVC of a 
stack of N = 9 coupled JJs with periodic BCs, at (3 = 0.6. 
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FIG. 14: (Color online) Breakpoint current Ibp and the jump 
current Ij (where transition to another CVC branch occurs) 
versus p. Inset shows the enlarged part where changes of the 
wave number of LPW happen. 



Inset (a) to this figure shows the enlarged BPR. In this 
inset we see that the chaotic part of BPR does not ap- 
pear. To stress the absence of the chaotic part in the 
BPR we also show the results for a system without noise 
(I n = 0) in Fig. [15] In this case the system behaves like 
a single J J without a BPR in CVC. The curves before 
the parametric resonance region coincide. To see the dy- 
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FIG. 15: (Color online) CVC of the stack with N=9 Josephson 
junctions at P — 0.6 and periodic BC. Inset (a) demonstrates 
the enlarged part of BPR. Inset (b) shows electric charge on 
layers versus time. 

namics of the electric charge in the S-layers, we show the 
time variation of charge on one layer of a stack with 9 



JJs with (3 = 0.6 and periodic BCs in the inset (b) to 
Fig. Q1)J We see absolutely regular oscillations. 

Finally, we look at the charge correlations in neighbor- 
ing S-layers in a system without chaotic behavior in the 
parametric resonance region. As we know such a study 
is a powerful tool for the investigation of the dynamics 
of coupled J Js^£. Results of calculations are presented in 
Fig. [TH]which demonstrates the charge-charge correlation 
function (as defined in Sec. Ill) versus the bias current 
for neighboring layers in the stack with 9 coupled JJs at 
/3 = 0.6. Here we are comparing the correlations within 
a specific stack, and so we use the usual notation for the 



correlation functions, i.e. C, 



1,1+1 



{1,1 + 1}. We see 
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FIG. 16: (Color online) (a) Charge-charge correlation func- 
tions versus bias current for the stack with 9 junctions for 
P = 0.6. Inset shows the charge distribution on S-layers at 
/ = 0.845. (b) Demonstration of correlation function pairing 
shown in (a). 

that the behavior of the correlation functions reflects the 
behavior of the CVC in the BPR shown in Fig. [131 an d 
corresponds to regular dynamics in the coupled JJs. The 
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correlations in the oscillations of the charge in S-layers 
persist up to the transition to another branch in CVC. 
The inset to Fig. ITST a) shows the charge distribution on 
S-layers at / = 0.845 and for the wave number k — An/%. 
In Fig. fTBT b) we see that the correlation functions form 
pairs corresponding to the standing LPW. This pairing is 
demonstrated in Fig. ll6f bL from which one can conclude 
that the 6th S-layer plays the role of the specific layer. In 
this case the charge on the 6th layer is non-zero, because 
the wavelength of the LPW is A = 4.5d (not shown in 
Fig. HSJb)). Thus, the above analysis is a proof of con- 
cept for the existence of parameters leading to a BPR 
without chaos. 



IX. CONCLUSIONS 

Physical properties of intrinsic Josephson junctions in 
high temperature superconductors continue to attract 
much attention due to the prospects they offer for su- 
perconducting electronics, in particular, as high power 
sources of coherent electromagnetic radiation in the THz 
region and high precision voltage standards. Applica- 
tions would depend strongly on the chaotic properties 
such systems. Observed radiation spectra from IJJs are 
very complex and are temperature dependent. Transition 
of any junction in the stack to the chaotic state would 
lead to a loss of synchronization and coherent emission 
from the stack. The dissipation parameter (McCumbcr 
parameter) is itself temperature dependent, and this has 
to be taken into account when predicting the transition 
to the chaotic state in real systems. 

Here we have studied chaotic features of coupled in- 
trinsic Josephson junctions which might impact on their 
possible applications. Different manifestations of the 
chaotic behavior were seen in the temporal dependence of 
the electric charge in the superconducting layers, phase- 
charge and charge-charge Lissajous diagrams, Lyapunov 
exponent and correlation functions. We demonstrated 
that the number of junctions in the stack and the bound- 
ary conditions do influence on the chaotic part in the 
breakpoint region. Experimental testing of the above fea- 
tures may therefore expedite progress in their potential 
applications. The chaotic part of BPR increases with the 
number of junctions in the stack, relative to the regular 



part of the BPR, at nonperiodic BCs. Transitions be- 
tween chaotic and regular behavior inside of the chaotic 
part of the breakpoint region were seen. It is anticipated 
that future developments in this direction of research may 
uncover new chaotic aspects in coupled JJs. 

Our analysis of the general features of the system has 
uncovered the conditions and parameters that produce a 
breakpoint region without chaos. We have shown that for 
high enough values of dissipation parameter j3 the chaotic 
part of breakpoint region does not exist. This is reminis- 
cent of a single overdamped junction which lacks chaotic 
dynamics.— We also arrived at the important prediction 
that the proximity effect should be expected to have a 
strong influence on the fine structure of the breakpoint 
region. 

In order to better understand the nature of the inter- 
mittency in coupled JJs the charge time dependencies 
in these regions should be considered in future work. It 
may also be of interest to calculate the full spectrum 
of Lyapunov exponents, rather than only the maximal 
exponent, in order to distinguish between chaotic and 
(possibly) hyperchaotic behavior in these systems. 

Presently there is renewed interest in phenomena re- 
lated to the switching from the outermost branch to 
inner branches and transitions between inner branches. 
Recently, for example, powerful coherent THz radiation 
was observed experimentally and was associated with a 
"bump" structure in the same part of the CVC where the 
BPR was found^iii The recently observed broadly tun- 
able sub-terahertz emission from internal branches of the 
CVCk 2 stresses the importance of further investigations 
of the chaos related to these inner branches. 
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